# RK spatial predictions of soil properties within individual sentinel sites
# e.g. Sindindi, Western Kenya
# Jiehua Chen & Markus Walsh, Dec. 2011

# set local working directory e.g.,
# setwd("~/Documents/LDSF/W. Kenya/Wet chem")

# load packages
require(nlme)
require(lattice)
require(proj4)
require(gstat)
require(geoR)

# to install the "rgdal" package on OSX use the following steps:
# download and install the GDAL framework from: http://www.kyngchaos.com/software/grass
# then in R specify:
# setRepositories(ind=1:2)
# install.packages("rgdal")
require(rgdal)

# Load soil profile coordinates and increment measurements
unzip("Sindindi_profiles.zip")
siinc <- read.table("Increment.csv", header=T, sep=",")
sipro <- read.table("Profile.csv", header=T, sep=",")

# Gridded 100 m res covariates (e.g SRTM Elevation, CTI, Slope, 2005 Landsat ETM+ bands 1-8, Landsat MrSID bands 2,4,7)
unzip("Sindindi_grids.zip")
grid.list <- c("cti.tif", "elev.tif", "etm1.tif", "etm2.tif", "etm3.tif", "etm4.tif", "etm5.tif", "etm61.tif", "etm62.tif", "etm7.tif", "pan.tif", "sid2.tif", "sid4.tif", "sid7.tif", "slope.tif")

# read "grid.list" into a "sigrid" data frame
sigrid <- readGDAL(grid.list[1])
names(sigrid)[1] <- sub(".tif", "", grid.list[1])
for (i in grid.list[-1]) {
	sigrid@data[sub(".tif", "", i[1])] <- readGDAL(paste(i))$band1
}

# project Lat/Lon profile coordinates in "sipro" to the UTM CRS of "sigrid"
sipro.utm <- as.data.frame(project(cbind(sipro$Lon, sipro$Lat), "+proj=utm +ellps=WGS84 +zone=36 +units=m +no_defs"))
colnames(sipro.utm) <- c("x", "y")
sipro <- cbind(sipro, sipro.utm)
coordinates(sipro.utm) <- ~x+y
proj4string(sipro.utm) <- CRS("+proj=utm +zone=36 +ellps=WGS84  +units=m +no_defs")

# overlay "sipro.utm" profile locations with "sigrid"
sipro.ov <- overlay(sigrid, sipro.utm)
sidat.ov <- cbind(sipro, sipro.ov@data)

# Topsoil/Subsoil data setup
sidat <- merge(sidat.ov, siinc, by="PID")
sitop <- sidat[sidat$TopSub=="topsoil",]
sisub <- sidat[sidat$TopSub=="subsoil",]

### aside, GLS models

# Generalized Least Squares (GLS) models
# GLS 1: Topsoil pH conditioned on toposequence covariates, assuming no spatial autocorrelation
# pHtop.gls <- gls(pH~I(cti-13)*I(elev-1321), data=sitop)

# GLS 2: GLS 1, updated to include spatial autocorrelation
# pHtop1.gls <- update(pHtop.gls, corr=corExp(form=~x+y, nugget=TRUE))
# summary(pHtop1.gls)

# GLS models 1 & 2 are nested so they can be compared them using the "anova" function
# anova(pHtop.gls, pHtop1.gls)

### end aside

## Regression Kriging (RK) models
# setup for "likfit" function in geoR
sitop.geo <- sidat[sidat$TopSub=="topsoil",]
sitop.geo$x <- sitop.geo$x/1000
sitop.geo$y <- sitop.geo$y/1000
sitop.geodata <- as.geodata(sitop.geo, coords.col=5:6, data.col=31, covar.col=7:8)

# REML 1: restricted max likelihood estimates in geoR
pHtop.reml <- likfit(sitop.geodata, ini=c(0.5, 0.5), trend=trend.spatial("cte", sitop.geodata, add=~I(cti-13)*I(elev-1321)), lik.met="REML")
summary(pHtop.reml)

# "krige" function setup in gstat
sitop <- sidat[sidat$TopSub=="topsoil",]
coordinates(sitop) <- ~x+y
proj4string(sitop) <- CRS("+proj=utm +zone=36 +ellps=WGS84 +units=m +no_defs")

# Linear model of pH conditioned on cti & elevation
pHtop.lm <- lm(pH~I(cti-13)*I(elev-1321), data=sitop)

# Variogram estimates estimate in gstat
pHtop.var <- variogram(residuals(pHtop.lm) ~ 1, sitop)
pHtop.vgm <- vgm(model="Exp", nugget=0.05, range=1000, kappa = 0.5)
pHtop.vgf <- fit.variogram(pHtop.var, model=pHtop.vgm)

# set variogram parameters equal to those from pHtop.reml
pHtop.vgf[1,2] <- pHtop.reml$nugget
pHtop.vgf[2,2] <- pHtop.reml$sigmasq
pHtop.vgf[2,3] <- pHtop.reml$phi*1000

# Variogram plot
plot1 <- plot(pHtop.var, pHtop.vgf, col="black")

# RK 1: Universal kriging estimates in gstat
pHtop.rk <- krige(pH~I(cti-13)*I(elev-1321), locations=sitop, newdata=sigrid, model=pHtop.vgf)

# Basic SP-plots
plot2 <- spplot(pHtop.rk["var1.pred"], scales=list(draw=TRUE), col.regions=rev(grey(seq(0,1,0.025))))
plot3 <- spplot(pHtop.rk["var1.var"], scales=list(draw=TRUE), col.regions=rev(grey(seq(0,1,0.025))))

# Export predictions to e.g. geotif
writeGDAL(pHtop.rk["var1.pred"], "pHtopEst.tif", type="Float32", driver="Gtiff")
writeGDAL(pHtop.rk["var1.var"], "pHtopVar.tif", type="Float32", driver="Gtiff")









